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TECHNICAL NOTE 2977 

TECHNIQUES FOR CALCULATING PARAMETERS OF NONLINEAR 
DTNAMIC SYSTEMS FROM RESPONSE DATA 
By Benjamin R. Briggs and Arthur L. Jones 


SUMMARY 


A study has "been made of the problem of determining nonlinear par- 
ameters of dynamic systems, such as aircraft or servomechanisms, whose 
measured responses are nonlinear. This problem cannot be solved in gen- 
eral but, for certain types of systems, adequate answers can be obtained. 
In the differential equations representing the motion of the latter 
systems, one or more coefficients will be functions of the nonlinear par- 
ameters. For the systems studied herein, the nonlinear parameters are 
functions of either the amplitude of the dependent variables or their 
time derivatives. The techniques developed for coefficient evaluation 
are, for the most part, extensions of two well-known methods of stabil- 
ity derivative evaluation for linear systems, namely, the derivative 
method and the method of inspection of transients for period and damping 
characteristics. Phase— plane methods such as those utilized in nonlinear 
mechanics are also included. 

It was found that the nonconstant coefficients could be determined 
satisfactorily for first— or second— order systems, using one or a combi- 
nation of the techniques developed. For higher— order systems the only 
practical means of obtaining results was found to be the derivative 
method. For some of the examples tried, however, this latter method 
gave poor results due to the occurrence of ill-conditioned algebraic 
equations. The n umb er of unknown coefficients, constant or nonconstant, 
to be determined also affected the accuracy of the results , the smaller 
the n umb er, of course, the better the chances for a satisfactory evalu- 
ation. 


INTRODUCTION 


Much attention has been given recently to the problem of determining 
stability derivatives of airplanes by means of flight and wind-tunnel 

dynamic testing. This problem of reducing response data to obtain stabil- 
ity par am eters is the inverse problem in dynamics, the direct problem 
being to find the response to an arbitrary forcing function when the 
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equations of motion anil the system parameters are known. For linear 
systems , methods for solving the inverse problem are summarized and 
evaluated in reference 1. More detailed descriptions of the individual 
methods can he found in references 2 through 6. 

In the investigation reported herein, the problem of determining 
nonlinear stability parameters from dynamic responses of nonlinear sys- 
tems was studied. Previous to performing the actual analytical opera- 
tions involved in this inverse type of analysis, it is desirable to 
determine the magnitude of the nonlinearities involved in order to class- 
ify the problem as either approximately linear or sufficiently nonlinear 
to necessitate the more complicated nonlinear approach. It is also nec- 
essary to make an assumption about the order and form of the differential 
equation of the system being investigated. The response of a system to 
sinusoidal inputs can be analyzed to find the magnitude of the nonline- 
arity involved, but there is no criterion for the identification of the 
order and form of a differential equation to represent the system. The 
latter must be determined more or less intuitively. 

There is a variety of effects which may cause nonlinear behavior 
in a dynamic system. These effects appear in several different ways in 
the equations of motion of the system. The cause of one common non- 
linear term is an element in the system whose response is dependent on 
the amplitude of a dependent variable. For example, it is possible for 
an airplane to have a nonlinear variation of pitching moment with angle 
of attack which would require a nonlinear term Cm(o-) to appear in the 
equations of motion in place of the usual linear term Cn^a. The tech- 
niques presented in this report are concerned essentially with systems 
having this type of nonlinearity. 

This investigation was mainly concerned with aircraft responses. 

The techn i ques presented, however, are shown to be applicable to responses 
of other mechanisms such as an autopilot servo system. The degree of non- 
linearity in the responses for high-speed aircraft of unorthodox configu- 
ration, such as a guided missile, is already a matter of concern to design 
engineers. (See ref. 7.) Furthermore, the automatic control systems 
employed in aircraft usually contain, within the servo system itself, 
nonlinearities due to factors such as amplifier saturation and velocity 
limitation. 

For the most part, the methods developed herein are extensions of 
certain linear techniques. In reference 1 the several methods for anar- 
lyzing linear responses are placed in two categories. The first cate- 
gory includes methods wherein the response data are substituted directly 
into the differential equations, the general forms of which are assumed. 
Such methods are ca ll ed "equations— of-motion methods." Examples of 
techniques in this category are the "derivative method" given in refer- 
ences 2 and 3> and the "matrix method" of reference 4. The second 
category involves methods wherein analytical curves are fitted to the 
dynamic responses. These techniques are known as "response curve— fitting 
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methods” and the methods described in references 5 and 6 are examples of 
this type. The same categorical divisions seem suitable for the nonlinear 
techniques derived from these linear ones. In addition, a technique of 
analysis based on phase— plane methods, which does not fit precisely into 
either of the above two categories, is presented for use in determining 
nonlinear damping functions associated with certain classes of self- 
sustained oscillations. 

Experience has shown that the results of the analysis of first— and 
second— order linear systems are more accurate than those of systems of 
higher order. It is not at all surprising, therefore, that this condi- 
tion also exists in the non li near situation. Fortunately, in many cases 
the dynamic behavior of systems can be adequately described by second- 
order differential equations whether linear or nonlinear. Emphasis has 
therefore been placed on techniques suitable for second— order systems. 

The derivative method, however, is applicable in principle in cases where 
the differential equations involved are of any order whatever, the pri- 
mary limitation being the accuracy in determining the higher-order deriv- 
atives of the response either by measurement or computation. In order to 
show that this limitation is not always serious, a successful demonstra- 
tion of the technique of analysis for a nonlinear system of fourth order 
with unknown parameter is presented in the illustrative examples. 


NOTATION 


j ma 


rate of change of pitching-moment coefficient with angle 
dCm 


of attach. 


5a 


C> 


za 


' z 5 




C m (a) 
Cm* (a) 


rate of change of normal— force coefficient with angle of 
attack , 


5c z 


5a 

rate of change of normal— force coefficient with elevator 

5C„ 

deflection angle, — - 

56 

rate of change of pitching-moment coefficient with pitching 

5c m 

velocity, — — 

5q 

nonlinear pitching-moment coefficient 

nonlinear pitching-moment parameter (See eq. (B15).) 
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Sa 


nonlinear undamped natural— frequency parameter (See eq. (Bl6).) 


P(x),F(i) 

G(x) 

H(t) 

I y 

Iz 

I(v e ) 

K 

S 
I % 


nonlinear damping coefficients, arbitrary second— order system 

nonlinear restoring— force coefficient, arbitrary second-order 
system 

forcing function, arbitrary second— order system 
moment of inertia about y axis 
moment of inertia about z axis 

nonlinear amplifier characteristic, autopilot servo 

I v pS 

nondimensional moment— of— inertia parameter, -4 $ — 


yawing moment 

rate of change of yawing moment with angle of yaw. 


Srr 

St 




rate of change of yawing moment with angle of sideslip. 


SN 


P 

PaAa 

Pf 

E 

S 

v o 

b 

S’ 


period of oscillation 

slope of nonlinear amplifier characteristic, autopilot servo 

gain parameter, autopilot servo 

radius of phase— plane limit cycles of Tan der Pol and Eayleigh 
equations, small values of p 

wing area 

free-stream velocity 

damping coefficient, second— order linear differential 
equation, 2^^ 

local wing chord 

1 r 3 —z> 

mean aerodynamic chord, — c^ dy 

no nli near damping term, Becond— order nonlinear differential 
equation 


f (x,x) 
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k 

kf»*m 

m 

q. 

t 

■fcl 


restoring— force coefficient, second— order linear differential 
equation, o^ 2 

parameters, autopilot servo 
mass 

pitching velocity 
time 

V 


■fcl/2 


Vf,Ti,V e 




y 

a 

0 

5 

e 


time to damp to half amplitude 

output, input, and error voltages of autopilot servo 

arbitrary dependent variables 

spauwise distance from plane of symmetry 

angle of attack 

angle of sideslip 

elevator deflection angle 

parameter, Whittaker’s smoothing formula 


P 

T 

T m> T D 

* 

A, A 2 , . . . 

03 

oun 

i 


density of air 
aerodynamic time, - - - -- 

psv 0 , 

parameters, autopilot servo 
angle of yaw 

parameter, Yan der Pol and Bayleigh equations 
first, second, etc., differences 
angular frequency of oscillation, radians/sec 
undamped natural frequency 
damping ratio 
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Subscript 


s smoothed, "by Whittaker’s method. 


SCOPE OF THE AHALYSIS 


The dynamic system under consideration, with one exception, will be 
assumed to be describable by a second— order differential equation such 
as 


x + f(x,±) + g(x) = H(t) 


( 1 ) 


where f(x,x) is either of the form F(x)x, or F(x)x, or by one or two 
first-order equations of the form 


x = f x (x,y,t) 
y = f 2 (x,y,t) 


( 2 ) 


In equations (l) and (2), x and y are arbitrary variables and the dots 
over x and y indicate differentiation with respect to time t. It 
will also be assumed that any nonlinear coefficients or parameters are 
functions only of the dependent variable or its first derivative, and 
not of the time. 

The exceptional case included will be a fourth-order system. This 
example is presented to demonstrate the application of the derivative 
method to nonlinear systems of higher order. The dynamic system in 
this example is an autopilot servo, and the differential equation 
describing its open— loop behavior is assumed to be of the form 


( 3 ) 


where Vf and v e represent output and error voltages and the coefficients 
Ao,Ax,A 2 , and A s are constants. The function l(ve) is nonlinear. 

In order to facilitate the actual numerical work involved in the 
illustrative examples presented in the appendixes, a Beeves Electronic 
Analogue Computer and IBM digital computing equipment were used. For the 
purpose of obtaining harmonic analyses of certain responses a harmonic 
analyzer was also employed. 
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EETERMDIAIIOIT OF DEffiEE OF NONLIHEABITY 


The degree of nonlinearity of a dynamic system, as determined "by 
examination of response data, is measured by the deviation of the 
responses from those of a purely linear system. It can he considered 
as an indication of the magnitude of the physical nonl i nearity present 
in a given system. For practical purposes it would he desirable to have 
some number, or figure of merit, which would be a quantitative measure 
of the intensity or degree of the no nli nearity. For example > a complex 
dynamic system may have a response that to all appearances is linear. 
However, if an analysis of the response data is to be attempted, it 
might save time and work if tests could be made to ascertain whether 
the data warranted linear or nonlinear treatment. All physical systems 
are essentially nonlinear to some extent, but if the degree of nonline- 
arity is small it is possible that linear analyses would suffice. 

There are certain response characteristics which are typical of 
nonlinear systems and which may be recognized without resorting to a 
figure of merit. One example is a system which possesses a self- 
sustained oscillation. As another example, the output of the system 
to two sinusoidal inputs of different frequency may possess not only the 
two impressed frequencies, but also multiples, and sums and differences 
of multiples of the impressed frequencies. These nonlinear phenomena 
and many more are discussed in detail in references 8 and 9. These 
effects are readily recognized in dynamic responses and they cannot in 
general be described by linear differential equations. When they are 
observed, then, there is no question but that nonlinear analyses must 
be used. 

On the other hand, there are cases where the effect of a nonline- 
arity in a dynamic system is not strongly evident upon cursory examina- 
tion of a response. A transient may have appearances of a linear 
response, for example, until it is observed that the period of the oscil- 
lation changes as the oscillation damps , which could be the effect of a 
nonlinear restoring force. The time to damp to half amplitude may be 
observed to vary for transient oscillations about different trim posi- 
tions, which is evidence of a nonlinearity in the damping. 

The presence or absence of nonlinearities in a dynamic system can 
usually be established by the examination of the ratio of output to 
input amplitude for various sinusoidal input magnitudes at a given 
frequency. If the amplitude ratio or the phase angle shows a substantial 
variation, the system is effectively non! inear. Such effects are most 
emphatically noted on frequency— response plots. These plots are com- 
monly presented with the logarithm of the amplitude ratio and the phase 
angle as ordinates and the logarithm of the frequency as the abscissa. 
(See, e.g., ref. 10.) For a purely linear system these curves do not 
vary with the amplitude of the input. In nonlinear systems, however, 
the resonance frequency may shift and peak amplitudes may increase or 
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decrease as the input amplitude varies. In the accompanying sketch 
logarithmic plots for a nonlinear system for three different amplitudes 



Log frequency 



Log frequency 


of a sinusoidal input are shewn. 

Any one of these curves has the 
appearance of a linear response hut 
since the characteristics of the 
system change as the amplitude of 
the input changes it must he con- 
cluded that the system is nonlinear. 

The deviations from linearity 
may he such that for engineering 
purposes, linear analyses could he 
used with sufficient accuracy. The 
simple techniques discussed in the 
preceding paragraphs, however, do 
not provide any figure of merit for 
so Judging. If limits are set for 
tolerable variation of such quanti- 
ties as time to damp to half ampli- 
tude, resonance frequency, peak 
amplitude, and phase shift for the 
range of inputs under consideration, 
then, it may he possible to estab- 
lish whether or not linear analyses 
can he used by examining transient 
responses or frequency— response data. 


As previously mentioned, however, it would he desirable to have a 
number or figure of merit indicative of the degree or intensity of non- 
linearity for a given system. It is possible to make a quantitative 
estimate of this sort, if steady— state responses to sinusoidal inputs 
are obtained. The figure of merit, called the distortion factor, 1 is 
defined as 100 times the square root of the sum of the squares of the 
amplitudes of all harmonics present in the response beyond the fund. a— 
mental, divided by the amplitude of the fundamental. For perfect sinus- 
oidal response data obtained from a purely linear system, the distortion 
factor is zero. Data obtained from sinusoidal inputs to nonlinear sys- 
tems, however, contain various amounts of higher harmonics, and the 
distortion factor is a quantitative measure of the amplitudes of these 
harmonics generated by the nonlinearities. 

In Appendix A a demonstration of the calculation of the distortion 
factor is presented. This quantitative test of nonlinearity was applied 
to data used subsequently in the examples illustrating the application 
of data reduction methods for nonlinear systems . The results indicated 
that a distortion factor in excess of about 5 percent would probably 
warrant application of nonlinear analyses. 

1 The distortion factor is commonly used in electrical and communica- 
tions engineering to measure wave distortion. 
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REDUCTION OF DATA TO NONLINEAR PARAMETERS 
Equations— of-Motion Method 


In reference 3 the derivative method, the only equations— of— 
motion technique that will he applied herein to nonlinear problems, is 
described for application to linear systems. This technique consists 
of obtaining the dynamic response of the system to a prescribed input, 
obtaining the proper number of response derivatives, and then for sev- 
eral instants of time substituting the corresponding values of input, 
response, and response derivatives into the differential equations of 
the system. In this manner a set of linear algebraic equations is 
developed wherein the variables are the coefficients of the differential 
equations. The coefficients are then calculated by solution of this set 
of algebraic equations either simultaneously or by least-squares methods. 
For linear dynamic systems a single response provides sufficient data 
for this calculation. A major difficulty in using this method of data 
reduction is in obtaining accurate derivatives of the output when they 
cannot be measured experimentally. 


The derivative method can be 
applied to reduction of nonlinear 
data in the following way. Sup- 
pose that several oscillations of 
the response of the system to a 
given input are recorded. For 
the sake of simplifying the fol- 
lowing discussion, let the given 
response be a free transient and 
assume that the system under con- 
sideration is of second order, 
describable by an equation of the 
general form of equation (l), 
where f(x,x) is of the form 
F(x)&. Obtain the first two deriv- 
atives of the response, then plot 
response and derivatives to the 
same time scale, as in the accom- 
panying sketch. At a value of x, 
say xo, draw a line parallel to 
the time axis on the x time his- 
tory. At each instant where this 
line cuts the response curve read 
the values of x and x. In this way 
several sets of values of x and x, 
all corresponding to the same value 
of x, namely x Q , are obtained. 

It should be pointed out here that 
there must be at least as many sets 
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of values of x and x as there are unknown coefficients in the dif- 
ferential equations of the system. 

Assume that the equation of motion of the system is of the form 


+ F(x)x + G(x)x = 0 (4) 

Also assume that for a given value of x, damping and. restoring force 
coefficients have unique values. These unique values will he taken as 
the numbers which are found hy solving hy least squares the set of 
equations which are obtained by substituting the sets of values of ± 
and. x corresponding to Xq, as determined by the graphical procedure 
of the previous paragraph, into equation (4). If this procedure is 
repeated for several more values of x, say X 3 _, x 2 , etc., then the 
numbers F(x 0 ), F(xx), etc., and G(x 0 ), G(x 2 .), etc., may be plotted 
against the corresponding values of x, and the nonlinear coefficients 
will then be represented graphically as functions of x. 

Consider now the application of the derivative method in the case 
where only one coefficient, suspected to be nonconstant, is unknown. 
From a single response of this system the unknown coefficient can be 
calculated by the derivative method. The coefficient will be valid 
within the range of values of displacement as given by the experimental 
response curve. To compute the coefficient, merely substitute the 
response and its derivatives into the differential equation and solve 
directly for the unknown parameter as a function of time. These results 
are then plotted against displacement to give the nonlinear variation 
of the parameter with displacement. 

Consider next the case where more than one coefficient is unknown. 
If it is assumed that all but one of these unknown coefficients are 
constant, there is a reasonable chance that the constant coefficients 
can be found by analyzing a small— amplitude response, using one of the 
linear techniques previously mentioned. However, if more than one 
coefficient must be solved for using the derivative method, there is a 
great possibility that the data used will result in ill-conditioned 
equations for the nonlinear systems. The apparent cause of this predi- 
cament is' that all the values of the derivatives of the dependent vari- 
able must be read at the one value of the dependent variable for one 
set of simultaneous equations. As a result, a l 1 the derivatives tend 
to have' nearly the same value in each equation. In order to avoid this 
difficulty it is desirable to use sets of data for which the particular 
value of the dependent variable being used occurs in widely separated 
parts of the response excursion. Such data could be obtained from 
sinusoidal inputs of different magnitudes or from a modulated sinusoidal 
input. A modulated input, wherein either the amplitude or frequency is 
modulated, has the advantage that the entire excursion range desired for 
the dependent variable could be covered in one input. The separate 
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sinusoidal inputs, however, were found to yield data that resulted in 
the less ill-conditioned simultaneous equations. 

As a general rule in applying the derivative method it is well to 
determine as many of the constant coefficients as possible, ‘beforehand, 
by other means. When the equations of motion of a system are derived 
there is one for each degree of freedom. It is preferable to use the 
several equations in the form in which they are derived, rather than to 
eliminate variables and obtain a single equation of high order. The 
necessity for taking derivatives of the output data will thereby be 
kept to a minimum. 

In Appendix B several numerical examples are presented in order 
to demonstrate the derivative method for one or two unknown parameters. 
In the first of these examples the numerical techniques for smoothing 
and differentiating response data are used and discussed in some detail. 
Examples 1 and 2 show that for a single unknown parameter very good 
results can be obtained, even for a system of fourth order. For the 
case of more than one unknown, however, demonstrated in the third 
example, the results were of poor accuracy due to the occurrence of ill- 
conditioned equations in the computations. For second— order systems 
with more than one unknown, better results may be obtained by rising the 
techniques presented in the next section. 


Response Curve-Fitting Methods 


Curve— fitting methods for response-data reduction are based upon 
the assumption that the responses can in some manner be approximated, 
or fitted, by analytical curves whose characteristics are known. One 
of the most elementary of the linear curve— fitting techniques is to 
examine transient responses to determine the period and damping charac- 
teristics, from which a sum of exponentials which fits the given tran- 
sient response can be written. A technique based upon this single 
linear transient-response analysis method, applicable to second— order 
nonlinear data, will be presented here. 2 

To carry out this technique of data reduction, responses of the 
second— order system are obtained to step inputs of varying amplitude. 
Now, the assumption is made that small oscillations about each of the 
steady-state, or trim, -displacements are linear, and a period and damp- 
ing analysis is made of each response to determine period P, and time 
to damp to half amplitude tj./ 2 . If the oscillation is too strongly 
damped, of course, this calculation is not feasible. Second— order 
linear differential equations are commonly written in the form 

£ The application of this well-known linear analysis method to nonlinear 
problems was suggested by Mr. Stewart M. Crandall of the Ames Aero- 
nautical Laboratory of the NACA. 
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x + 2 ^DhX + a)n 2 x = H(t) 


(5) 


(see, e.g. , ref. 10) where £ is the damping ratio and £% is the 
natural undamped frequency. It is also well known that when H(t) is 
a step or pulse input, P and tx/ 2 are related to C and in the 
following manner: 


■fci/2 


0.693 




P = 


2rt 

CD 


' 


( 6 ) 


a> = % 




The values of P and t x / 2 for each of the step input responses are 
then substituted into equations (6) to determine the coefficients 
2£o> q anrl cdq 2 , valid, under the present assumptions, at the particular 
trim value of the displacement. 



The values of c% 2 and 2^c% 
are then plotted against the trim 
displacement to which they correspond, 
thus determining their variation with 
displacement. The accompanying 
sketches illustrate this graphical 
procedure. 

To demonstrate this technique, 
a numerical example is presented in 
Appendix B. The determination of 
both the damping coefficient and a 
restoring-force term of a second- 
order system by this method is found 
to be relatively easy and accurate. 

A somewhat simpler technique 
which can be used only for finding 
restoring— force nonlinearities, based 
upon inspection of transient responses, 
can be developed from the steady-state 
method given in reference 3. For a 
step input to a second— order system, 
the response assumes a constant steady- 
state value, where the derivatives of 
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the response and tlie input are all zero. Thus, in the steady state, 
the differential equation is reduced to a relationship between the step 
input and the restoring force. By using this technique on the first 
system considered in Appendix B (see eq. (A3)), for example, (^(a) can 
he determined from the following relationship: 


C m (a 0 ) = i__ (2.285a D - 4.6205 o ) 

66.181 


where ao is the steady— state value of a, resulting from a step input 
of magnitude S 0 . In this manner, the nonlinear parameter can he deter- 
mined as a function of steady— state values of the response, as was done 
previously in a more elaborate fashion. The advantage of this technique 
is that aperiodic responses can he analyzed, hut a disadvantage is that 
the damping c ann ot he calculated directly. 

Another curve— fitting method based on linear equivalence, so 
obvious that it is usually the first to he suggested, was tried. In 
this method, the linear damping or restoring-force equivalents obtained 
from excursions centered around zero displacement and covering various 
amplitudes of the response would he used to construct the nonlinear 
function. Unfortunately, however, no basis was found for translating 
these equivalent linear slopes into the correct nonlinear function. 

The simple technique of fairing through the end points of the resulting 
family of linear slopes did not provide a satisfactory result. 


Phase^Plane Methods, Self-Sustained Oscillations 


In this section, an inverse process will he presented for the cal- 
culation of the nonlinear damping term in the differential equation of 
a second— order system which exhibits self-sustained oscillations. The 
snaking oscillations of an aircraft are self— sustained, for example, 
and it is possible that the nature of the nonlinear damping of such 
oscillations could he determined by phase— plane methods. It should he 
mentioned, however, that unless the damping corresponds fairly closely 
to the nonlinear damping of either the Van der Pol or Bayleigh classi- 
cal equation of nonlinear mechanics , the phase— plane method presented 
here may not he useful. 

For any oscillation having a periodic variation of the displace- 
ment with time, the time rate of change of displacement is likewise 
periodic. Thus, the phase-plane plot, displacement versus velocity, 
for such an oscillation is a closed curve. Moreover, if the closed 
curve is unique, that is, if the trajectory of the motion in the phase 
plane eventually reaches the closed curve regardless of initial condi- 
tions, the oscillation has a limit cycle. The limit cycle is a 
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characteristic of self— sustained oscillations and is markedly different 
from the behavior of forced or undamped linear systems -which can also 
exhibit closed trajectories in the phase plane. In these latter cases, 
however, the closed trajectories vary in size and shape for differing 
initial conditions. 

It is known that a large class of second— order systems capable of 
self-sustained oscillations can be described by either a Van der Pol 
equation or a Rayleigh equation. In the former equation the damping 
term is of the form 


f(y,y) = -(a - By 2 ) f (7) 

and in the latter the damping term is 

f(y,30 = -(Cy - Dy 3 ) (8) 



The nature of the damping coefficient 
of Van der Pol’s equation is shown 
in the first of the accompanying 
sketches, and the nature of the damp- 
ing term of Rayleigh’s equation is 
shown in the second sketch. In gen- 
eral, the restoring force need not be 
linear for a limit cycle to exist 
(see ref. 11). For the present anal- 
ysis, it is assumed that the observed 
self— sustained oscillation can be 
described by a Van der Pol or a 
Rayleigh equation, and that the restor- 
ing forces are linear. 

When, in equation (l), the restor- 
ing force is linear and the damping 
force is of the form of equation ( 7 ) 
the result is Van der Pol’s equation 

y - (A - By 2 ) y + o^y = 0 ( 9 ) 

If the transformations 
ti = (Pat 1 



A 


x 


( 10 ) 
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are made, then equation (9) becomes 

- n(l - x 2 ) i + 2 = 0 (11) 

dti dt i 

Equation (30) is the form in which Van der Pol’s equation is .usu ally 
written. Similarly, if Eayleigh’s equation, 


y _ (Cy _ Dy 3 ) + o^ 2 y - 0 


( 12 ) 


is transformed hy means of 

t x = o^t 



The result is 


& 2 x _ „ [ dx_ _ 1 / dx V 
dt x 2 litx 3 \cLti ) 


+ x = 0 


(13) 


m 


This is the canonical form of Eayleigh’s equation. Equations (11) and 
(l4) are related by simple transformations, which will be shown next. 

/ \ dz 

In equation (11 J, let x = , and the result is 

dti 



d 2 z dz 
dt^ dti 


= 0 


(15) 


Integrate equation (15) and the result is 


d 2 z 

dti 2 


dz 1 
_dtl “ 3 



- n 


+ z = const. 


(16) 
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If tlie change of variable z — const. = x is made, equation (l6) 

■becomes identical to equation (14). Thus, the Yan der Pol equation has 
"been transformed into a Rayleigh equation. Conversely, the Eayleigh 
equation, equation (l4) , can he transformed into a Yan der Pol equation 

hy the change of variable z = . The result of this change of vari— 

dti 

able is differentiated once, which puts it into the form of equation (ll). 
These manipulations show that the solution of Yan der Pol’s equation is 
the derivative of the solution of Eayleigh’ s equation. 

From theoretical considerations (see ref. 8, pp. 197—199)/ it can 
he shown that for small amounts of damping (i.e., small values of p 

or — the phase— plane limit cycle of Yan der Pol’s equation is almost 

a circle of radius 



and the limit cycle of Eayleigh’ s equation "with small values of p 

n 

or -— is nearly a circle of radius 

* (i8) 

When Yan der Pol’s and Rayleigh’s equations are written in the standard 
forms, equations (ll) and (l4), then the phase plots are almost circles 
with radii R = 2 for small values of p. 

In reference 12 a series expansion valid for all values of p, is 
given for the solution of a Yan der Pol equation in the form of equa- 
tion (ll). From this series it can he shown that the maximum value of 
displacement does not change appreciably for Increasing values of p, 
hut remains nearly 2. By means of the simple relationship between 
Yan der Pol’s and Eayleigh’s equations, and the above mentioned series 

solution, it can he shown that the variable -^2- in equation (l4), 

dti 

Eayleigh’s equation, changes hut slightly as p increases. This prop- 
erty is utilized in the inverse process, which will he described next. 

The essence of the inverse technique is to compare the phase-plane 
plot of the experimental data to previously prepared phase plots of the 
standard forms of Yan der Pol’s or Eayleigh’s equation. Such plots were 
obtained hy solving equations (ll) and (l4) on the analogue computer for 
several values of p. (See figs. 1 and 2.) It will he observed in 
figures 1 that the maxi mum displacement does remain sensibly constant 
as p increases, and in figures 2 that the maxi trmm value of the 
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derivatives of the displacement does not change appreciably as p 
increases. Note that figures 1 are quite different in appearance from 
figures 2, particularly at large values of p. 

Suppose that a dynamic system is assumed descrihahle by a 
Tan der Pol equation. The frequency of oscillation of the response 
can be observed from the experimental responses, and it is approxi- 
mately equal to c%. Thus a value of is obtained. Next the 

phase plot of x versus — is prepared. The maximum value of x is 

a>n 

noted, and it is set equal to E (eq.(17)). This fixes a value for 
the ratio S.. The next step is to plot the transformed phase portrait, 

the variables of which are 



This new phase plot is compared to the previously prepared phase plots 
of equation (ll), figures 1, from which an approximate value of p is 
determined. Now, using the relation 




A_ 


( 20 ) 


a value of A can be computed. With this value, and equation (17), B 
can be calculated. Thus the coefficients of the Van der Pol type of 
equation which approximately fits the experimental data are found. 


If the given data are assumed to be responses of a Bayleigh equa- 
tion, the technique for finding the coefficients of the damping term is 
similar to that of the foregoing discussion. A phase plot, y versus 

Y 

is made and the maximum value of — is observed. This value is then 

r 

related to the ratio through equation (18). The transformed phase 
portrait, the variables being % 



dx _ / 3Iton 2 jr_ 
dti J C 


( 21 ) 
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is drawn nest, and it may "be compared to figures 2 to ascertain an 
approximate value of p. Now 


IX = ± (22) 

so C is dete rmine d . Substituting C into equation (l8) will give a 
value for D. In this manner C and D are computed for the system 
which, is describable "by a Rayleigh equation. 

An example showing application of the foregoing technique is pre- 
sented in Appendix B. The responses used in this demonstration were 
obtained from wind-tunnel tests of an aircraft model mounted on a bal- 
ance that permitted free oscillation in yaw and roll. Some difficulty 
is encountered in fixing a precise value of p in this problem, due to 
the fact that fluctuating conditions in the wind tunnel resulted in the 
responses not being exactly periodic. The manner in which a compromise 
value of p is chosen is given in the above-mentioned appendix. 


DATA AND METHODS ENQUIRED EC® D1TERMIMTI0N 
OF NONLINEAR PARAMETERS 


The procedure for reducing nonlinear response data to values of 
the parameters of the physical system being studied begins with knowl- 
edge or assumptions concerning the order and form of the differential 
equations describing the system. Such information can usually be 
obtained by careful consideration of the mechanical or electrical com- 
ponents of the system. Once suitable assumptions regarding the differ- 
ential equations have been made, the nature of the various parameters 
can be investigated by methods such as those discussed and referred to 
herein. 


First-Order Systems 


When the mass of a one-degree— of— freedom mechanical system is 
negligible with respect to the damping- andrestoring— force coefficients, 
the system can be approximately represented by a first-order differen- 
tial equation. Such a situation miay also occur in an electrical sys- 
tem having very small inductance. These are known as degenerate sys- 
tems. Thiq parameters of degenerate systems can be determined by the 
derivative method, using one or more dynamic responses. 
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Second-Order Systems 


The derivative method can he applied in the case of second— order 
systems with any or all of the parameters unknown. The data recommended 
for the calculation of several parameters hy this method are responses 
to sinusoidal inputs of various magnitudes, or modulated sinusoidal 
input responses. It was found, in the calculations given in a previous 
section, that the difficulty of ill-conditioned algebraic equations was 
not greatly reduced when responses to modulated inputs were used, how- 
ever. The most satisfactory kind of data appeared to he several 
responses to sinusoidal inputs having greatly differing amplitudes. 

When the restoring— force and damping— force coefficients are the 
unknowns, and if either or both are nonlinear, then the curve— fitting 
technique wherein responses to several step inputs are used can he 
applied with much more assurance of success than can the derivative 
method. It is of interest to note that this curve— fitting technique 
is independent of parameters in the forcing function, which may or may 
not he known. 

If any one parameter, such as the restoring-force coefficient, 
the damping coefficient, or a coefficient in the forcing function is 
unknown, then the derivative method may he applied. A single dynamic 
response, such as a sinusoidal, or pulse or step input response is all 
the data that are needed. In the special case of nonlinear damping, 
which is partly dissipative and partly regenerative, limit cycles can 
occur, and this sort of nonlinearity is best studied hy means of phase- 
plane techniques. One or more cycles of the self— sustained vibration 
of the response are the required data. 


Systems of Higher Order Than Second 


In all cases of this type the derivative method seems to he the 
only practical approach. The number of responses required, preferably 
to sinusoidal inputs of different amplitudes, will depend entirely on 
the number of unknown parameters. 


CONCLUDING REMARKS 


Methods have been presented for reducing response data from certain 
nonlinear dynamic systems to coefficients or system parameters. The non- 
constant coefficients were assumed to he functions either of the depend- 
ent variable or of its time derivative. The methods could he separated 
into three categories, namely, the equations-of-motion methods, the 
response curve— fitting methods, and the phase-plane methods. 
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Of the first category, the derivative method, which was applicable 
to all the cases investigated, was found to he most suitable when only 
one unknown nonconstant coefficient was to be determined. As the number 
of unknown coefficients, constant or nonconstant, increased, the accu- 
racy of the derivative method decreased due to the ill-conditioned 
nature of the s imult aneous equations to be solved. Furthermore, as the 
order of the differential equation representing the system increased, 
the accuracy of the derivative method tended to decrease due to the 
difficulty of measuring or computing higher— order derivatives of the 
dependent variable. 

The response curve— fitting method employed was the most satisfactory 
for second-order systems having more than one unknown coefficient. If 
the damping force is known to be partly degenerative and partly regen- 
erative, such that self-sustained oscillations are observed, the phase- 
plane method was found to be a suitable method of determining the form 
of the nonlinear damping in certain types of second— order systems. 


Ames Aeronautical Laboratory 

National Advisory Committee for Aeronautics 
Moffett Field, Calif. , Apr. 14, 1953 
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APPENDIX A 

CALCULATION OF THE DISTORTION FACTOR 


The purpose of this appendix is to demonstrate numerical techniques 
for determining the distortion factor of the responses of a dynamic 
system to sinusoidal inputs. The dynamic system under consideration 
here could he characteristic of a missile, and only the longitudinal 
motion of the missile will he studied. The equations of motion of this 
system, which are used in the numerical examples of Appendix B also, 
are presented here. 

The longitudinal equations of motion of the misBile are assumed in 
the form 


2m — 2Tq — Cz^ct — — Cz^6 
2 Kt q — Cin(a) — Cm^q = Um§S 


> 


(Al) 


For a particular situation the coefficients in equation (Al) were calcu- 
lated to he 


r 


°“q 

° z 5 

Cm- 


= 2.188 sec 

= 0.001578 
=-0.0284 sec 

= 0.030 

= 0.070 

= -5.32 


The numerical yalues of the above 
coefficients are based upon a 
being measured in radians, q in 
radians per second, and 5 in 
degrees. The function C m (a) is 
nonlinear and of the form shown in 
the sketch. Equations (Al) can now 
he written 


a-q+1. 2157a + 0.006855 = 0 
q+1.8795q - 66.1810^(0,) - 4.63275 =0 


(A2) 
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If q is eliminated "between the two equations (A2), the result is 

a + 3.0952& + 2.2850a - 66.l8lC m (a) = 4.6198& - O.OO 685 S (AB) 


The data for the present computation were obtained by solving 
equation (A3) on the analogue computer with 5 = S 0 sin art where to 
is 2ic radians per second, and where values of 5 0 used were 0.5 > 1°, 

2°, 3°, and 3.5°. The 5 term was not included. The function Cm(a) 
as determined from wind-tunnel tests is shown in figure 3 . The responses, 
to the sinusoidal inputs are presented in figure 4. 

By use of a wave analyzer, the amplitudes of the first 10 harmonics 
of a steady-state cycle of each of the above responses were found. The 
amplitude of the first harmonic was adjusted to the value 100 so that 
the computation of the distortion factor would be simplified. The 
results are summarized in the chart that follows: 


Harmonic 

Relative amplitudes of harmonics 




an 

&o = 3.f 

1 

100.00 

100.00 

100.00 

100.00 

100.00 

2 

8.95 

13.02 

16.00 

13.90 

12.90 

3 

3.61 

2.08 

4.85 

6.23 

7.90 

4 

1.14 

2.00 

1.09 

1.68 

1.68 

5 

.37 

.26 

.20 

1.42 

1.75 

6 

;54 

.67 

.20 

.32 

.44 

7 

.42 

.12 

.25 

.54 

.64 

8 

.15 

.24 

.25 

.21 

.55 

9 

.11 

.35 

.33 

.42 

.50 

10 

.10 

.33 

.16 

.47 

.27 

Distortion 

factor 

9.8 

13.3 

16.7 

15.4 

15.3 


As a check on the accuracy of the wave analyzer two exact sine waves 
were analyzed. The distortion factors for the two cases were 2.0 
and 1.5. 

By 'use of the same steady— state responses (fig. 4) the harmonic 
content of each was calculated by a digital integration process with 
the aid of the IBM computing equipment. The relative amplitudes of 
the first five harmonics of each response cycle are presented in the 
following chart: 
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Harmonic 

Relative amplitudes 

of harmonics 

Wmm 

IBS 

0 

OJ 

n 

0 

<0 

5o = 3° 

8 0 = 3.5° 

1 

100.00 

100.00 

100.00 

100.00 

100.00 

2 

9.85 

15.68 

19.10 

17.25 

17.21 

3 

.87 

1.30 

5.51 

5.62 

7.14 

h 

1.16 

1.14 

.81 

1.90 

.90 

5 

.29 

.65 

.51 

.48 

1.17 

Distortion 

factor 

10.0 

15.8 

19.8 

18.3 

18.6 


An exact sine wave was analyzed in this manner and the percentage dis- 
tortion was nearly zero. 

The distortion factors were smaller when the wave analyzer was used 
than when the digital analysis was employed. There is probably no 
particular significance in this trend. It may be said, however, that 
there is much less chance for error in the digital process since the 
data to be used are read directly from the analogue computer records. 

For the wave analyzer it was necessary to cut a template in the shape 
of the wave to be analyzed, and there existed the probability of intro- 
ducing inaccuracies at this step. 

There are many methods for computing harmonic coefficients by hand. 
(See, for instance, refs. 13 and 14. ) In the digital method used in 
the IBM computation, 2h points were taken over each of the steady-state 
response cycles. Then parabolas were fitted to overlapping triples of 
these points, thus giving 12 parabolas for each response. The integrals 
which define the Fourier coefficients were then integrated exactly, 
using these parabolic approximations', to the function being analyzed, 
and the 12 integrals were summed to obtain each Fourier coefficient. 

This technique is quite accurate, but is probably too lengthy for hard 
computation. From the results of the previous analyses, it is apparent 
that the amplitude of the first harmonic divided by the amplitude of the 
fundamental closely approximates the value of the distortion factor. 

Thus any of the harmonic analyses which may be performed quickly by 
hand on desk calculating machines or slide rules probably give a good 
indication of the distortion factor, if enough points are used to 
assure a reasonably accurate determination of the first two harmonics 
beyond the fundamental. 
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APPENDIX B 

ILIESTRAirTE EXAMPLES 

Numerical Examples of the Derivative Method 


Inasmuch as the success of the derivative method depends to a great 
extent upon the accuracy of the derivatives of the response, which are 
usually computed rather than measured, the numerical techniques employed 
herein for the data differentiation will he discussed briefly. The 
differentiations were performed by means of finite difference formulas, 
which are described in detail in reference 13. Since experimental data 
contain errors and irregularities, it is necessary to smooth the data 
mathematically before a numerical differentiation can be performed with 
accuracy. In the present work this smoothing has been carried out by 
means of a technique which is presented in reference 13 , pages 303 — 316 , 
and in references 15 and 1 6 . Examples of the application of this method 
of smoothing, to be called herein the Whittaker method, will be fo und in 
these references. 

In order to illustrate the effect of the Whittaker smoothing 
method, data which have been smoothed will be compared to the unsmoothed 
data in example 1. It will be shown that a primary effect of smoothing 
is that the third differences of the data are much smaller and more 
regular than the third differences of the unsmoothed data. The differ- 
entiation formulas therefore converge quite rapidly by using only the 
first three differences of the smoothed data. 

Since the smoothing process is quite tedious, it has been per- 
formed on IBM computing equipment. The differentiations were all 
carried out by means of desk— type calculating machines, however. 

Example 1. determination of nonlinear restoring— force parameter , 
all other parameters known, of a second-order system .— With the aid of 
the analogue computer, equation (A3) was solved using as the input 
5 = sin 2rtt and with Cm(a) given by the cubic 


C m (a) = -1.2a- 4.0a 2 - 90a 3 (Bl) 

The response a is shown in figure 5. This response is now regarded 
as the given data, and (^(a). is to be calculated. All of the other 
parameters of equation (A3) will be assumed known. 

If several cycles of the response were available in the steady- 
state range , then an average of corresponding ordinates for several 
such cycles would undoubtedly reduce some of the irregularities in the 
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data. In the present example, however, there were not sufficient data 
for this, hut the first five points were averaged with the first five 
points of the following cycle so that there would he more assurance 
that the numerically computed derivatives of the assumed cycle of data 
would appear as derivatives of a truly periodic response. 

Ordinates of a were read from figure 5 at intervals of 0.02 
second, starting at t = 1.84 seconds and ending at t = 2.92 seconds. 
These data are recorded in column 1 of table I. The first and last 
five values were then paired and averaged in the manner mentioned in 
the preceding paragraphs. These averaged end values are given in 
column 2 of table I. With these new end ordinates the data were assumed 
to he a representative steady-state cycle of a. The first three dif- 
ferences of the response a were taken, and in column 3 of table I the 
third differences were recorded, and for purposes of comparison later 
on, the derivative da/dt was taken and recorded in column 4. 

Whittaker’s smoothing method was then applied to the a data, 
using 0.05 as the value of the parameter e in the smoothing formula. 
The smoothed values are denoted by the symbol ag. In table I and 
columns 5 , 6, and 7> these values are tabulated, along with their third 
differences and the averaged end values of Os . Note that the third 
differences of the smoothed data, designated by the symbol A 3 as, are 
much smaller than those of the unsmoothed data, while the actual devia- 
tion of smoothed from unsmoothed values is very slight. 

The as data were next differentiated and the result was denoted 
dcxg/dt. These data are tabulated in table I, column 8. After smooth- 
ing da s /dt, the result being designated by the symbol (dOg/dt) g , and 
after pairing and averaging the first and last five values of the 
smoothed data, a second derivative d(dag/dt) s /dt was taken. For com- 
parative purposes, the derivative of the unsmoothed data da s /dt was 
also computed. In columns 9> 10, 11, and 12 of table I the data per- 
taining to the second derivative of a are tabulated. 

In figure 6 the derivatives of a and Og are plotted, and in 
figure 7 the derivatives of da s /dt and (dOs/dt) s are likewise shown. 
These comparisons serve to show the need for, and effect of, the smooth- 
ing of a set of data prior to performing a numerical differentiation. 

The derivatives that have been computed will now be used in con- 
nection with equation (A3) to demonstrate the derivative method of data 
reduction. 

The forcing function of equation (A3) is 


H(t) = 4.6198 5 - 0.00685 5 


(B2) 
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which, is, in this example, 


H(t) = 4.6198 sin 2«t -0. Oh-31 cos 2*t (B3) 

Solving equation (A3) for (^(a) gives 


C m (a) = — i [a + 3.0952 a + 2.2850 a - H(t) ] (B4) 

66.181 

In equation (b 4) the data of columns 5, 10, and. 12 of table I were used 
for a, a, and a, respectively. The forcing function, given by equa- 
tion (B2), was tabulated in column 13 of table I. By use of these 
data Cm(c0 was first computed as a function of time. Since each value 
so found corresponds to a value of a, a plot of C m (t) versus a(t) gave 
the desired function, C m (a). In column l4 of table I Cm(a) is given. 

The function Cm(a) is plotted against a in figure 8. If the 
results were regarded as mathematically exact, then it might be con- 
cluded that Cm (a) is double valued. It is known, however, that C m (a) 
is single valued, so this apparent hysteresis is, in fact, due to sys- 
tematic distortion injected into the calculation in the smoothing and 
differentiating processes. For a more direct comparison of this result 
with the known C m (a) (eq. (Bl)), a cubic was fitted by least-squares 
methods to the computed function. This cubic is 


C m (a) = -1.234 a - 4.004 a 2 - 40.456 a 3 (B5) 

The large discrepancy between the coefficients of a 3 in equations (B5) 
and (Bl) is 'undoubtedly caused by the strong effects of a 3 which are 
only felt at large values of a, and the data of this example do not 
extend sufficiently far from zero to afford a more accurate value- for 
this coefficient. Note in figure 8 that the function C m (a) of equa- 
tion (B5) is quite accurate except at large values of a. 

As a further example in computing a single unknown parameter by 
the derivative method, the techniques described in the first example 
will be applied to transient responses of equations (A2). The data 
were obtained by solving equations (A2) on an analogue computer, using 
for Cm(a-) the function given in figure 3. The oscillation was initi- 
ated by an initial displacement in a of 8° (0.1396 radians). The 
responses were thus free oscillations, and the forcing term 6 was not 
present in the computations to determine (^(a).. 
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In figure 9 are shown portions of the responses q and a of the 
above system. The unknown function C m (a) can be calculated from the 
second of equations (A 2 ), so only one derivative of q was required, 
and none of a. The responses in q and a were read at intervals of 
0.02 second for a 1-second range starting at t = 0, and these results 
are given in table II. The q data were smoothed by Whittaker's method, 
with e => 0.1, and the first derivative was taken. These results are 
also presented in table II. In figure 10 the derivative dq g /dt is 
plotted, and it is apparent that the value of e of 0.1 is not too 
large. The accuracy with which the unsmoothed a and q can be read 
from the analogue computer traces in this case justifies the use of a 
value larger than was used in the first example. 

The function Cjafa) was found from the equation 


0m(a) - sriar (i + 1 - 8795 q) 


(B6) 


The results of this calculation are given in table II and are plotted 
in figure 11. 

Example 2. determination of amplifier characteristic curve in auto- 
pilot servo system, all other parameters known.— In reference 17 the 
dynamic behavior of an autopilot servo system was simulated on the ana- 
logue computer. An amplifier in the servo system was known to possess 
a saturation type of nonlinear characteristic. In order to account for 
certain time lags in the dynamic system it was found necessary to intro- 
duce an exponential lag operator (see ref. 10) into the equations of 
motion. A three— term approximation was made to the lag operator, which 
had the effect of raising the differential equation from second to fourth 
order. By use of this fourth^-order differential equation and experimen- 
tal responses obtained in bench tests of the autopilot, the problem in 
this example was to calculate the nonlinear amplifier characteristic 
when all the other parameters of the system are' known. 


After the approximation to the lag operator has been made, the open- 
loop equation is 


T D 2T m d \ 


Vm + 




2Pfkfkjn dt 


, 3 2 

d Vf Tp+T m 4 Yf 

+ + 


dvf 


Bfkfkjn dt 3 Pfkfkm dt 2 Pfkfkm dt 


I(v e ) (B7) 


where vf is the output voltage and v e the error voltage. The error 
and output are related to the input by the relation v e = vj_ — vf , 
where Vj, is the input voltage. The parameters xp, T m , Pf, k f , 
and k^ are all constants, and 
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l(y e ) = ’PoAaJe (B8) 

characteristic . 


(B9) 


When these values are used, equation (B7) becomes 

d 4 Vp d 3 V-p d 2 Vp d vp , 

o.oooon + 0.00263 + 0.315 + 5.17 i(t 6 ) = o (bio) 


where P a Aa is the nonconstant slope of the amplifier 
In a certain case the parameters in equation (l4) were 

P f = 0.24 

kf = 12.8 

k m = O.O63 [ 

r m = 0.052 

t D = 0.009 J 


The response of the servo system to an input voltage 
Vi = 1.56 sin 2itt was obtained. In figure 12, vp, Vf, and v e are 
plotted. The function Vf> , does not appear at a glance to be the 
response of a nonlinear system. By use of the digital harmonic analysis 
technique described in Appendix A, the distortion factor was computed 
and found to be 7.68« 

After smoothing Vf> by Whittaker’s method with e = 0.25, the 
first derivative of Vf> was found, using the difference formulas previ- 
ously cited. Three more derivatives were calculated, smoothing each 
one before the next derivative was taken. These four derivatives are 
shown in figure 13 . It is apparent that even with the large magnitudes 
attained by the fourth derivative, its contribution in equation (17) 
was very small. The four derivatives of v^ were substituted into 
equation (17). Then l(v e ) was computed in the same manner as was 
Cm (a,) in the previous two examples. The results of this calculation 
and all the other data for this example are given in table III. In 
figure 14, l(v e ) is plotted as calculated here and also as found by 
static tests of the autopilot. The false hysteresis effect is present, 
but the actual saturation trend of the characteristic is preserved in 
the function as computed by the derivative method. 

Example 3, determination of two unknown parameters by the deriva— 
tive method.— The extension of the derivative method to the analysis of 
n onU nflflr rasporiRe data for two unknown coefficients will now be demon- 
strated with numerical examples. The data were obtained by solving 
equation (A3) on the analogue computer, where Cm (a) is again given by 
figure 3, and the 6 term is neglected. For the purpose of the follow- 
ing computations, equation (A3) was put into the form 
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a + ba + kx = 0 o 5 


(Bll) 


where b, k, and C Q are constants. Equation (Bll) is therefore linear, 
hut it was used to describe the behavior of the. system under considera- 
tion only over small ranges of values of a. The coefficient b was 
assumed known from other considerations, and k and 0 o were calculated 
as constants for several values of a, using the technique previously 
described. The computed values of k and C Q could then be related to 
the values of a at which they are determined in order to ascertain 
their variation as a is varied. In this particular example the func- 
tion C m (a) was calculated in addition to k, using the relation 

C m (a) = a (B12) 


Equation (B12) was obtained from the restoring-moment terms of equa- 
tion (A3), and k is the number computed above for the value of a at 
which Cn^a) was to be determined. 

A set of data was obtained for the phase-modulated sinusoidal input 


B = 1.75 sin (2rtt + « sin 0.3«t) 


(B13) 


As was mentioned in a previous section, this sort of response data has 
the advantage that one run of the dynamic system is sufficient. The 
data so obtained were smoothed and- differentiated. Then the response, 
its first two derivatives, and the forcing function were plotted to the 
same time scale in. the manner described previously. At values of a 
of 0.02, 0.06, and 0.08 radians, lines were drawn through the response 
curve parallel to the time axis. At each instant where a line cuts 
the a curve, the corresponding values of the two derivatives and the 
forcing function were read. For each value of a, then, the quantities 
so read were substituted into equation (Bll). The resulting sets of 
equations were solved by least squares for k and C Q . These calcula- 
tions are summarized in the following table, where k is also related 
to C m (a) by equation (B12): 


a 

Co 

k 

Cm (a) from 
equation (B12) 

Cnx(O') from 
figure 3 

0.02 

4.861 

111'. 42 

-0.033 

-0.030 

.06 

6.082 

161.51 

-.144 

-.120 

.08 

5.446 

150.31 

-.179 

-.180 
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The correspondence between the calculated values of 0 ^( 0 , ) and the 
values taken from figure 1 is indicated in the above chart, anfl the 
true value of C Q is 4.62. 

As a second example, an amplitude-modulated sine wave was used 
for 8; that is. 


6=2 sin 2itt + sin 2.2rtt (Bl4) 

Again b was assumed known, and calculations similar to these above 
were made to determine values of k and C Q for several values of a. 
These results are given in the following table: 


a 

Co 

k 

C m (a) from 
equation (B12) 

(^(a) from 
figure 3 

0.025 

2.669 

49.25 

- 0.018 

-0.032 

.050 

3.806 

101.22 

-.075 

-.095 

-.075 

5.154 

105.9^ 

.118 

.110 


Besponses of this same system to sinusoidal inputs of various 
amplitudes were obtained, and a calculation similar to the ones dis- 
cussed above was made. The values of C Q were more accurate in the 
case of the responses to sinusoidal inputs, but Cm (a) was of about 
the same order of accuracy as before. 

The results of the above calculations are not entirely satisfactory. 
The difficulty of iU-canditi oning of the least-squares normal equations 
was present, although modulated inputs were tried. 

From these examples it is apparent that the accuracy of the deri- 
vative method of nonlinear data reduction may be poor for the case 
of two unknown parameters, due to the seemingly unavoidable ill- 
conditioned equations that must be solved. For systems of order 
higher than second, however, this method may yield results which can 
then be further refined. 
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Numerical Example of the Curve-Fitting Technique 


In order to demonstrate the simple period and damping technique of 
analysis of no nli near response data, the following example is presented. 
The data were obtained by solving equation (A3), with the & term 
omitted, on the analogue computer, using step inputs of magnitudes -0.5°, 
-1°, -2°, -3°, - 4 °, O. 25 0 , 0.5°, 0.75°, 1°, 2°, 3°, and 4 °. The func- 
tion C m (a) that was used is that given in figure 3. With the assump- 
tion that each response is linear, ti/ 2 and P were computed for the 
responses. It was found that ti/ 2 is nearly constant, and the 
values of 2 £<%, the damping coefficient, are tabulated here. The aver- 
age value for 2£wn is 2.91., which differs by about 5 percent from the 
correct constant value, 3.095. This indicates that the technique gives 
reasonably accurate results. 


It was observed that the per- 
iod P of the oscillations varies 
as the amplitude of the input varies. 
Thus 0^2 varies. In order to make 
a comparison with Cm (a) of the var- 
iation of oi^ 2 let 


0^0 co = ■ — 66 . 181 


8Cm* 

~£T a 


(B 15 ) 


The numbers 



are calculated 


8a 

and then plotted against the corre- 
sponding trim values of a. This curve, 
given in figure 15 , was then integrated 
in an approximate manner and the result was 
bol CjnMoO. This function is related to 


5 o 

25 % 

4.0 

3 .U 

3.0 

3.00 

2.0 

2.91 

1.0 

2.89 

.75 

2.85 

.5 

2.85 

.25 

2.78 

-.5 

3.15 

-1.0 

3.00 

-2.0 

2.89 

-3.0 

2.96 

— 4.0 

2.89 


designated by the sym- 
Cm(a) by the equation 


-66.181 C m * (a) = 2.285 a - 66.181 Cm(a) (Bl6) 


The function C m (cO, as calculated from equation (Bl6), is plotted in 
figure 16. The Cm(a-) of figure 3 is also given in figure l6 for com- 
parison. Considering the basic simplicity of this calculation, the 
results are very good. 
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Numerical Example of the Phase-Plane Method 


As an example of the application of the phase-plane technique for 
determining the nonlinear damping involved in a self— sustained oscilla- 
tion, a numerical problem will be solved. In wind-tunnel tests per- 
formed upon an aircraft model, the model was observed to initiate a 
self— sustained lateral oscillation about its mounting as the speed of 
the wind in the tunnel reached a certain level. Several cycles of this 
oscillation were recorded and read at intervals of 0.005 second, and 
these data are tabulated in table IV. The above data were smoothed by 
Whittaxer * s method with e = 0.25, and then differentiated. The results 
of this calculation are also presented in table TV. 

The differential equation of the dynamical system was assumed to 
be of the form 


I z *+ f(iM0+ (% + Np)* = 0 (B17) 

The constants I z and (N^ + Np) were known to be 

I z = 0.9^ in— lb— sec 2 /radian 
N^ + Np = 7^2.5 in-lb /radian 

and f(^,^) was to be calculated. If the notation of the previous dis- 
cussion is adhered to, then 


a>n = 



28. 05 radians/sec 


The phase portrait, versus $/(£a, was drawn next. (See fig. 17. ) 
Note that the phase curve is not a closed limit cycle. This is due to 
variations in the data, caused undoubtedly by fluctuating conditions 
in the wind tunnel. From the general trend of this phase portrait, 
namely t he predominantly negative slopes in the first and third quad- 
rants, flufl also fram^physical considerations that the damping should be 
a function only of i}r, it was concluded that a Bayleigh equation can be 
used to describe the motion. The damping term of equation (B17) was 
therefore assumed to have the form 
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= _(c if- Df 3 ) (B18) 

Iz 

From figure 17 the average maximum value of */cD n vas found to "be 
approximately 

K = 0.391(57.3) = 22.4a 

The conversion factor 57.3 is required, since f is given in degrees 
and On has the dimensions of radians per second, and K must he 
dimensionless. Then, using equation (18), it was found that 


J = 36.^0 


(B19) 


By means of equations (21), the new transformed phase portrait was con- 
structed. The variables were 


x = 5.10* 


dx llr 

= 5.10 — 

dt x <°n 


(B20) 


This phase plot (fig. 18) was compared to figure 2, and it was found 
that the phase curve of p =0.6 fits the experimental phase portrait. 
This curve is superimposed upon figure 18 for comparison. From the 
value of p of 0.6 and the relation p = C/cDq the value of C was 
found to he 16.82, and then from equation (B19)> D was calculated. 

It was 612. With these two numbers, equation (B17) can he written 


0.944t- (15.9*- 578* ) + 742.5* = 0 


(B21) 
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TAB IE I.- CALCULATIONS FOE EXAMPLE 1 
[Determination of (^(a) from, response to sinusoidal input] 
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TABLE II, _ CALCULATION'S ECB EXAMPLE 1 
[Determination of (^(a) from transient response data ] 
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TABLE III.- CALCULATIONS FOE NONLINEAR SERVO EXAMPLE 
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Figure 1.— Limit cycles of Van der Pol's equation for various values of 


S' 

H 




■25 -2j 0 -15 -t.0 -05 0 05 10 15 2.0 25 

Displacement, x 

(c) n = 0.2 


Figure 1.— 


25 -2 JO -15 -1.0-05 0 05 t.0 

Displacement, x 

(oL) n m 0.4 

rued. 


15 2.0 25 



Velocity,— 



(e) (J. » 0.6 


(f) H = 0.8 


Figure 1.— Continued.. 
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Figure 1.— Continued. 
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Figure 1.— Continued. 


vn 



■25 -2D -15 -LO -a5 0 0.5 LO 15 25 25 

Displacement, x 

(k) ij. = 4.8 


-25 -2D - 15 -tO -05 O 05 LO L5 20 25 

Disnlacement. x 
, / ” 

(l) p. n 8.0 


Figure 1.— Concluded.. 
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Figure 2.— Limit cycles of Bay! 
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Figure 2.— Continued. 
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Figure 2.— Continued. 
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Figure 2.— Continued. 







Figure 3." Nonlinear variation of pitching-moment coefficient with angle of attack from w ind — 

tunnel tests of a missile. 
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Figure 10.— Numerically computed, derivative of the smoothed transient pitching velocity of 

a missile. 
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Figure 11.— Comparison of true nonlinear pitching-moment coefficient curves with, curve calculated 

by the derivative method. 
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Figure 13.— Tlie first four derivatives with respect to time of the output 

voltage of a fourthr-order servo system. 
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Figure 17.— Phase— plane plot of the self— sustained oscillation of an air- 
craft model in a wind tunnel. 
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Figure 18.— Comparison of the phase— plane plot of the self— sustained 
oscillation of an aircraft model and the limit cycle of Rayleigh’s 
equation with p =0.6. 
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